This week we’ll cover a few packages for visualizing our data:

library(tidyverse)
library(timetk)
library(plotly)
library(readxl)
library(janitor)
library(lubridate)

knitr::opts_chunk$set(message = FALSE, warning = FALSE, comment = NA)

Let’s load our HPI dataset, wrangle it, and pivot it into long format.

url <- "https://www.fhfa.gov/DataTools/Downloads/Documents/HPI/HPI_PO_monthly_hist.xls"

destfile <- "hpi.xls"

curl::curl_download(url, destfile)

hpi <- read_excel(destfile, skip = 3)

hpi_wrangled <- hpi %>% 
  clean_names() %>% 
  slice(-1) %>%   # remove empty row
  rename(date = month) %>%
  select(date, ends_with("_sa")) %>%  # only keep seasonally adj. data
  # separate() from tidyr package to split date into separate columns for day/month/year 
  separate(date, into = c("year", "month", "day"), sep = '-', convert = TRUE, remove = FALSE) %>% 
  # unite() from tidyr to join columns 
  unite(yr_mon, year, month, sep = "/", remove = FALSE) %>% 
  mutate_if(is.numeric, round, digits = 2) %>%  # round all numeric columns to 2 digits
  # add labels with case_when()
  mutate(season = case_when(between(month, 3, 5) ~ "spring",
                            between(month, 6, 8) ~ "summer",
                            between(month, 9, 11) ~ "fall",
                            # between(month, 12, 2) ~ "winter" == won't work because there's no numbers between 12 and 2
                            TRUE ~ "winter")) %>% 
  select(date, yr_mon, year, month, day, season, everything()) %>%  # reorder columns
  # arrange() from dplyr to sort rows
  arrange(date)

hpi_tidy <- 
  hpi_wrangled %>% 
  select(date, contains("north"), contains("south")) %>% 
  # pivot_longer makes data long, or tidy
  pivot_longer(-date, names_to = "division", values_to = "hpi") %>% 
  group_by(division)

hpi_tidy

ggplot2

And now let’s create a basic ggplot.

ggplot(hpi_tidy, aes(x = date, y = hpi, color = division)) +
  geom_line()

The ggplot2 package uses some unique syntax (the “grammar of graphics”) that allows us to create highly customizable static graphics. This grammar can be a bit hard to grasp, so don’t worry if it takes a while to “click”.

I think about ggplot like this:

Let’s dissect the code above to understand each step.

Wait, why a + instead of %>%? You can think of ggplots in layers. At each step, we’re adding a new layer, like we’re painting on a canvas. This is different than the pipe, which is for passing an object along to a new function.

Faceting

Another thing we can “add” to our ggplots is a faceting layer. Facets divide a plot into subplots based on one of our variables. For example:

ggplot(hpi_tidy, aes(x = date, y = hpi, color = division)) +
  geom_line() +
  facet_wrap(~division)

Scatterplot

There are many kinds of geoms we can add. For example, a scatterplot uses geom_point().

hpi_pct <- 
  hpi_tidy %>% 
  mutate(pct_change = (hpi / lag(hpi)) - 1,
         pct_change_12_mons = (hpi / lag(hpi, 12)) - 1) %>%
  na.omit()

hpi_pct %>% 
ggplot(aes(x = pct_change_12_mons, y = pct_change, color = division)) +
  geom_point() + #alpha = .5
  facet_wrap(~division, ncol = 4) +
  labs(x = "% change (1 mo.)", y = "% change (1 yr.)") +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        plot.title = element_text(hjust = 0.5))

Scatterplot with trend line

We can also add multiple geoms to a ggplot object, e.g. adding a trend line to our scatterplots:

hpi_pct %>% 
ggplot(aes(x = pct_change_12_mons, y = pct_change, color = division)) +
  geom_point() + #alpha = .5
  geom_smooth(method = "lm",  se = TRUE, color = "purple") +
  facet_wrap(~division, ncol = 4) +
  labs(x = "% change (1 mo.)", y = "% change (1 yr.)") +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        plot.title = element_text(hjust = 0.5))

Histogram

There are lots of ways to customize the look and feel of your plots:

hpi_pct %>% 
  ggplot(aes(x = pct_change)) +
  geom_histogram(fill = "darkblue", color = "darkred", bins = 50) +
  facet_wrap(~division) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  theme_minimal()

Density

hpi_pct %>% 
  ggplot(aes(x = pct_change)) +
  geom_density(fill = "darkblue", color = "darkred", bins = 50) +
  facet_wrap(~division) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  theme_minimal()

Histogram and Density

hpi_pct %>% 
  ggplot(aes(x = pct_change)) +
  geom_histogram(fill = "darkblue", color = "darkred", bins = 20) +
  geom_density(color = "red") +
  facet_wrap(~division) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

Shading for recession

We can even combine datasets on the same plot. Let’s create a geom to shade the background of a plot during recession periods.

# First, create our recessions tibble with tribble()
recessions <- 
tribble(
  ~Peak, ~Trough,
  "1948-11-01", "1949-10-01",
  "1953-07-01", "1954-05-01",
  "1957-08-01", "1958-04-01",
  "1960-04-01", "1961-02-01",
  "1969-12-01", "1970-11-01",
  "1973-11-01", "1975-03-01",
  "1980-01-01", "1980-07-01",
  "1981-07-01", "1982-11-01",
  "1990-07-01", "1991-03-01",
  "2001-03-01", "2001-11-01",
  "2007-12-01", "2009-06-01",
  "2020-02-01", "2020-05-01"
  ) %>% 
  mutate(Peak = ymd(Peak),
         Trough = ymd(Trough))


recession_shade <- 
  geom_rect(data = recessions, 
            inherit.aes = F, 
            aes(xmin = Peak, 
                xmax = Trough, 
                ymin = -Inf, 
                ymax = +Inf), 
            fill = 'pink', 
            alpha = 0.5)


hpi_pct %>% 
  ggplot(aes(x = ymd(date), y = pct_change, color = division)) +
  recession_shade +
  geom_line() +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust=0)) +
  ylab("") +
  xlab("Percent change, monthly") +
  ggtitle("Housing Price Appreciation", 
          subtitle = "by US Census Division") +
  labs(caption = "data source: FHFA") +
  scale_x_date(limits = c(as.Date(min(hpi_pct$date)), as.Date(max(hpi_pct$date)))) 

timetk for time series data

The timetk package includes a bunch of functions that make working with time series data super easy. This includes functions for easily creating great looking plots of time series data:

hpi_pct %>% 
  ungroup() %>% 
  plot_time_series(date, pct_change, .color_var = division, .smooth = FALSE)

Anomaly diagnostics with timetk

timetk also includes functions for automatic anomaly detection:

hpi_pct %>% 
  filter(division == "south_atlantic_sa") %>% 
  ungroup() %>% 
  plot_anomaly_diagnostics(date, pct_change)

Plotly and ggplotly()

timetk’s interactive plots rely on plotly, a library for building interactive JavaScript visualizations. Plotly is supported in several different languages (including R) and has its own syntax.

Importantly, the plotly R package includes a function called ggplotly() that (you guessed it!) turns ggplots into plotly charts.

LS0tCnRpdGxlOiAiRGF0YSBWaXN1YWxpemF0aW9uIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpUaGlzIHdlZWsgd2UnbGwgY292ZXIgYSBmZXcgcGFja2FnZXMgZm9yIHZpc3VhbGl6aW5nIG91ciBkYXRhOgoKKiBbYGdncGxvdDJgXShodHRwczovL2dncGxvdDIudGlkeXZlcnNlLm9yZy8pLCBvdXIgbWFpbiBwYWNrYWdlIGZvciBzdGF0aWMgZGF0YSB2aXN1YWxpemF0aW9uIHdpdGggdGhlIHRpZHl2ZXJzZQoqIFtgdGltZXRrYF0oaHR0cHM6Ly9idXNpbmVzcy1zY2llbmNlLmdpdGh1Yi5pby90aW1ldGsvaW5kZXguaHRtbCksIGEgcGFja2FnZSBidWlsdCBieSBNYXR0IERhbmNobyBmb3Igd29ya2luZyB3aXRoIHRpbWUgc2VyaWVzIGRhdGEKKiBbYHBsb3RseWBdKGh0dHBzOi8vZ2l0aHViLmNvbS9yb3BlbnNjaS9wbG90bHkpLCBhbiBpbnRlcmFjdGl2ZSB2aXN1YWxpemF0aW9uIGJ1aWxkZXIgdGhhdCBhbHNvIGludGVncmF0ZXMgbmljZWx5IHdpdGggZ2dwbG90LgoKYGBge3Igc2V0dXAsIG1lc3NhZ2UgPSBGQUxTRSwgd2FybmluZyA9IEZBTFNFfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeSh0aW1ldGspCmxpYnJhcnkocGxvdGx5KQpsaWJyYXJ5KHJlYWR4bCkKbGlicmFyeShqYW5pdG9yKQpsaWJyYXJ5KGx1YnJpZGF0ZSkKCmtuaXRyOjpvcHRzX2NodW5rJHNldChtZXNzYWdlID0gRkFMU0UsIHdhcm5pbmcgPSBGQUxTRSwgY29tbWVudCA9IE5BKQpgYGAKCgpMZXQncyBsb2FkIG91ciBIUEkgZGF0YXNldCwgd3JhbmdsZSBpdCwgYW5kIHBpdm90IGl0IGludG8gbG9uZyBmb3JtYXQuCgpgYGB7cn0KaHBpIDwtIHJlYWRfZXhjZWwoIkhQSV9QT19tb250aGx5X2hpc3QueGxzIiwgc2tpcCA9IDMpCgpocGlfd3JhbmdsZWQgPC0gaHBpICU+JSAKICBjbGVhbl9uYW1lcygpICU+JSAKICBzbGljZSgtMSkgJT4lICAgIyByZW1vdmUgZW1wdHkgcm93CiAgcmVuYW1lKGRhdGUgPSBtb250aCkgJT4lCiAgc2VsZWN0KGRhdGUsIGVuZHNfd2l0aCgiX3NhIikpICU+JSAgIyBvbmx5IGtlZXAgc2Vhc29uYWxseSBhZGouIGRhdGEKICAjIHNlcGFyYXRlKCkgZnJvbSB0aWR5ciBwYWNrYWdlIHRvIHNwbGl0IGRhdGUgaW50byBzZXBhcmF0ZSBjb2x1bW5zIGZvciBkYXkvbW9udGgveWVhciAKICBzZXBhcmF0ZShkYXRlLCBpbnRvID0gYygieWVhciIsICJtb250aCIsICJkYXkiKSwgc2VwID0gJy0nLCBjb252ZXJ0ID0gVFJVRSwgcmVtb3ZlID0gRkFMU0UpICU+JSAKICAjIHVuaXRlKCkgZnJvbSB0aWR5ciB0byBqb2luIGNvbHVtbnMgCiAgdW5pdGUoeXJfbW9uLCB5ZWFyLCBtb250aCwgc2VwID0gIi8iLCByZW1vdmUgPSBGQUxTRSkgJT4lIAogIG11dGF0ZV9pZihpcy5udW1lcmljLCByb3VuZCwgZGlnaXRzID0gMikgJT4lICAjIHJvdW5kIGFsbCBudW1lcmljIGNvbHVtbnMgdG8gMiBkaWdpdHMKICAjIGFkZCBsYWJlbHMgd2l0aCBjYXNlX3doZW4oKQogIG11dGF0ZShzZWFzb24gPSBjYXNlX3doZW4oYmV0d2Vlbihtb250aCwgMywgNSkgfiAic3ByaW5nIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGJldHdlZW4obW9udGgsIDYsIDgpIH4gInN1bW1lciIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICBiZXR3ZWVuKG1vbnRoLCA5LCAxMSkgfiAiZmFsbCIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAjIGJldHdlZW4obW9udGgsIDEyLCAyKSB+ICJ3aW50ZXIiID09IHdvbid0IHdvcmsgYmVjYXVzZSB0aGVyZSdzIG5vIG51bWJlcnMgYmV0d2VlbiAxMiBhbmQgMgogICAgICAgICAgICAgICAgICAgICAgICAgICAgVFJVRSB+ICJ3aW50ZXIiKSkgJT4lIAogIHNlbGVjdChkYXRlLCB5cl9tb24sIHllYXIsIG1vbnRoLCBkYXksIHNlYXNvbiwgZXZlcnl0aGluZygpKSAlPiUgICMgcmVvcmRlciBjb2x1bW5zCiAgIyBhcnJhbmdlKCkgZnJvbSBkcGx5ciB0byBzb3J0IHJvd3MKICBhcnJhbmdlKGRhdGUpCgpocGlfdGlkeSA8LSAKICBocGlfd3JhbmdsZWQgJT4lIAogIHNlbGVjdChkYXRlLCBjb250YWlucygibm9ydGgiKSwgY29udGFpbnMoInNvdXRoIikpICU+JSAKICAjIHBpdm90X2xvbmdlciBtYWtlcyBkYXRhIGxvbmcsIG9yIHRpZHkKICBwaXZvdF9sb25nZXIoLWRhdGUsIG5hbWVzX3RvID0gImRpdmlzaW9uIiwgdmFsdWVzX3RvID0gImhwaSIpICU+JSAKICBncm91cF9ieShkaXZpc2lvbikKCmhwaV90aWR5CmBgYAoKIyBgZ2dwbG90MmAKCkFuZCBub3cgbGV0J3MgY3JlYXRlIGEgYmFzaWMgYGdncGxvdGAuCgpgYGB7cn0KZ2dwbG90KGhwaV90aWR5LCBhZXMoeCA9IGRhdGUsIHkgPSBocGksIGNvbG9yID0gZGl2aXNpb24pKSArCiAgZ2VvbV9saW5lKCkKYGBgCgpUaGUgYGdncGxvdDJgIHBhY2thZ2UgdXNlcyBzb21lIHVuaXF1ZSBzeW50YXggKHRoZSAiZ3JhbW1hciBvZiBncmFwaGljcyIpIHRoYXQgYWxsb3dzIHVzIHRvIGNyZWF0ZSBoaWdobHkgY3VzdG9taXphYmxlIHN0YXRpYyBncmFwaGljcy4gVGhpcyBncmFtbWFyIGNhbiBiZSBhIGJpdCBoYXJkIHRvIGdyYXNwLCBzbyBkb24ndCB3b3JyeSBpZiBpdCB0YWtlcyBhIHdoaWxlIHRvICJjbGljayIuCgpJIHRoaW5rIGFib3V0IGdncGxvdCBsaWtlIHRoaXM6CgoqIFdlIGFsd2F5cyBzdGFydCBieSBjYWxsaW5nIGBnZ3Bsb3QoKWAgdG8gY3JlYXRlIG91ciBwbG90IG9iamVjdC4gVGhpcyBpcyAodXN1YWxseSEpIHdoZXJlIHdlIHdpbGwgc3BlY2lmeSBvdXIgZGF0YSBhbmQgKmFlc3RoZXRpYyBtYXBwaW5ncyosIG9yIGhvdyBvdXIgdmFyaWFibGVzIHNob3VsZCBtYXAgb250byBmZWF0dXJlcyBvZiB0aGUgcGxvdCBsaWtlIGF4ZXMsIGNvbG9yZWQgZ3JvdXBpbmdzLCBldGMuCiogVGhlbiB3ZSBhZGQgKmdlb21zKiB0byBvdXIgcGxvdC4gVGhpcyBpcyB0aGUgc3RlcCB0aGF0IHdpbGwgYWN0dWFsbHkgZGlzcGxheSBvdXIgZGF0YSBvbiB0aGUgZ3JhcGguIFdlIHdpbGwgYWx3YXlzIGhhdmUgYXQgbGVhc3Qgb25lIG9mIHRoZXNlLCBhbmQgc29tZXRpbWVzIG11bHRpcGxlLgoqIEZpbmFsbHksIHdlIGNhbiBjaGFuZ2UgdGhlIGFwcGVhcmFuY2Ugb2Ygb3VyIHBsb3QgYnkgYWRkaW5nIHRoaW5ncyBsaWtlICp0aGVtZXMqIGFuZCBjaGFuZ2luZyB0aXRsZXMgYW5kIGNhcHRpb25zLgoKTGV0J3MgZGlzc2VjdCB0aGUgY29kZSBhYm92ZSB0byB1bmRlcnN0YW5kIGVhY2ggc3RlcC4KCldhaXQsIHdoeSBhIGArYCBpbnN0ZWFkIG9mIGAlPiVgPyBZb3UgY2FuIHRoaW5rIG9mIGdncGxvdHMgaW4gbGF5ZXJzLiBBdCBlYWNoIHN0ZXAsIHdlJ3JlIGFkZGluZyBhIG5ldyBsYXllciwgbGlrZSB3ZSdyZSBwYWludGluZyBvbiBhIGNhbnZhcy4gVGhpcyBpcyBkaWZmZXJlbnQgdGhhbiB0aGUgcGlwZSwgd2hpY2ggaXMgZm9yIHBhc3NpbmcgYW4gb2JqZWN0IGFsb25nIHRvIGEgbmV3IGZ1bmN0aW9uLgoKCiMjIEZhY2V0aW5nCgpBbm90aGVyIHRoaW5nIHdlIGNhbiAiYWRkIiB0byBvdXIgZ2dwbG90cyBpcyBhIGZhY2V0aW5nIGxheWVyLiBGYWNldHMgZGl2aWRlIGEgcGxvdCBpbnRvIHN1YnBsb3RzIGJhc2VkIG9uIG9uZSBvZiBvdXIgdmFyaWFibGVzLiBGb3IgZXhhbXBsZToKCmBgYHtyfQpnZ3Bsb3QoaHBpX3RpZHksIGFlcyh4ID0gZGF0ZSwgeSA9IGhwaSwgY29sb3IgPSBkaXZpc2lvbikpICsKICBnZW9tX2xpbmUoKSArCiAgZmFjZXRfd3JhcCh+ZGl2aXNpb24pCmBgYAoKCgojIyBTY2F0dGVycGxvdAoKVGhlcmUgYXJlIG1hbnkga2luZHMgb2YgZ2VvbXMgd2UgY2FuIGFkZC4gRm9yIGV4YW1wbGUsIGEgc2NhdHRlcnBsb3QgdXNlcyBgZ2VvbV9wb2ludCgpYC4KCmBgYHtyfQpocGlfcGN0IDwtIAogIGhwaV90aWR5ICU+JSAKICBtdXRhdGUocGN0X2NoYW5nZSA9IChocGkgLyBsYWcoaHBpKSkgLSAxLAogICAgICAgICBwY3RfY2hhbmdlXzEyX21vbnMgPSAoaHBpIC8gbGFnKGhwaSwgMTIpKSAtIDEpICU+JQogIG5hLm9taXQoKQoKaHBpX3BjdCAlPiUgCmdncGxvdChhZXMoeCA9IHBjdF9jaGFuZ2VfMTJfbW9ucywgeSA9IHBjdF9jaGFuZ2UsIGNvbG9yID0gZGl2aXNpb24pKSArCiAgZ2VvbV9wb2ludCgpICsgI2FscGhhID0gLjUKICBmYWNldF93cmFwKH5kaXZpc2lvbiwgbmNvbCA9IDQpICsKICBsYWJzKHggPSAiJSBjaGFuZ2UgKDEgbW8uKSIsIHkgPSAiJSBjaGFuZ2UgKDEgeXIuKSIpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIsCiAgICAgICAgYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA5MCwgaGp1c3QgPSAxLCB2anVzdCA9IDAuNSksCiAgICAgICAgcGxvdC50aXRsZSA9IGVsZW1lbnRfdGV4dChoanVzdCA9IDAuNSkpCmBgYAoKCiMjIFNjYXR0ZXJwbG90IHdpdGggdHJlbmQgbGluZQoKV2UgY2FuIGFsc28gYWRkIG11bHRpcGxlIGdlb21zIHRvIGEgYGdncGxvdGAgb2JqZWN0LCBlLmcuIGFkZGluZyBhIHRyZW5kIGxpbmUgdG8gb3VyIHNjYXR0ZXJwbG90czoKCmBgYHtyfQpocGlfcGN0ICU+JSAKZ2dwbG90KGFlcyh4ID0gcGN0X2NoYW5nZV8xMl9tb25zLCB5ID0gcGN0X2NoYW5nZSwgY29sb3IgPSBkaXZpc2lvbikpICsKICBnZW9tX3BvaW50KCkgKyAjYWxwaGEgPSAuNQogIGdlb21fc21vb3RoKG1ldGhvZCA9ICJsbSIsICBzZSA9IFRSVUUsIGNvbG9yID0gInB1cnBsZSIpICsKICBmYWNldF93cmFwKH5kaXZpc2lvbiwgbmNvbCA9IDQpICsKICBsYWJzKHggPSAiJSBjaGFuZ2UgKDEgbW8uKSIsIHkgPSAiJSBjaGFuZ2UgKDEgeXIuKSIpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIsCiAgICAgICAgYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA5MCwgaGp1c3QgPSAxLCB2anVzdCA9IDAuNSksCiAgICAgICAgcGxvdC50aXRsZSA9IGVsZW1lbnRfdGV4dChoanVzdCA9IDAuNSkpCmBgYAoKCgojIyBIaXN0b2dyYW0KClRoZXJlIGFyZSBsb3RzIG9mIHdheXMgdG8gY3VzdG9taXplIHRoZSBsb29rIGFuZCBmZWVsIG9mIHlvdXIgcGxvdHM6CgpgYGB7cn0KaHBpX3BjdCAlPiUgCiAgZ2dwbG90KGFlcyh4ID0gcGN0X2NoYW5nZSkpICsKICBnZW9tX2hpc3RvZ3JhbShmaWxsID0gImRhcmtibHVlIiwgY29sb3IgPSAiZGFya3JlZCIsIGJpbnMgPSA1MCkgKwogIGZhY2V0X3dyYXAofmRpdmlzaW9uKSArCiAgdGhlbWUoYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA5MCwgaGp1c3QgPSAxLCB2anVzdCA9IDAuNSkpICsKICB0aGVtZV9taW5pbWFsKCkKYGBgCgoKIyMgRGVuc2l0eQoKYGBge3J9CmhwaV9wY3QgJT4lIAogIGdncGxvdChhZXMoeCA9IHBjdF9jaGFuZ2UpKSArCiAgZ2VvbV9kZW5zaXR5KGZpbGwgPSAiZGFya2JsdWUiLCBjb2xvciA9ICJkYXJrcmVkIiwgYmlucyA9IDUwKSArCiAgZmFjZXRfd3JhcCh+ZGl2aXNpb24pICsKICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDkwLCBoanVzdCA9IDEsIHZqdXN0ID0gMC41KSkgKwogIHRoZW1lX21pbmltYWwoKQpgYGAKCgojIyBIaXN0b2dyYW0gYW5kIERlbnNpdHkKCmBgYHtyfQpocGlfcGN0ICU+JSAKICBnZ3Bsb3QoYWVzKHggPSBwY3RfY2hhbmdlKSkgKwogIGdlb21faGlzdG9ncmFtKGZpbGwgPSAiZGFya2JsdWUiLCBjb2xvciA9ICJkYXJrcmVkIiwgYmlucyA9IDIwKSArCiAgZ2VvbV9kZW5zaXR5KGNvbG9yID0gInJlZCIpICsKICBmYWNldF93cmFwKH5kaXZpc2lvbikgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gOTAsIGhqdXN0ID0gMSwgdmp1c3QgPSAwLjUpKQpgYGAKCgojIyBTaGFkaW5nIGZvciByZWNlc3Npb24KCldlIGNhbiBldmVuIGNvbWJpbmUgZGF0YXNldHMgb24gdGhlIHNhbWUgcGxvdC4gTGV0J3MgY3JlYXRlIGEgZ2VvbSB0byBzaGFkZSB0aGUgYmFja2dyb3VuZCBvZiBhIHBsb3QgZHVyaW5nIHJlY2Vzc2lvbiBwZXJpb2RzLgoKYGBge3J9CiMgRmlyc3QsIGNyZWF0ZSBvdXIgcmVjZXNzaW9ucyB0aWJibGUgd2l0aCB0cmliYmxlKCkKcmVjZXNzaW9ucyA8LSAKdHJpYmJsZSgKICB+UGVhaywgflRyb3VnaCwKICAiMTk0OC0xMS0wMSIsICIxOTQ5LTEwLTAxIiwKICAiMTk1My0wNy0wMSIsICIxOTU0LTA1LTAxIiwKICAiMTk1Ny0wOC0wMSIsICIxOTU4LTA0LTAxIiwKICAiMTk2MC0wNC0wMSIsICIxOTYxLTAyLTAxIiwKICAiMTk2OS0xMi0wMSIsICIxOTcwLTExLTAxIiwKICAiMTk3My0xMS0wMSIsICIxOTc1LTAzLTAxIiwKICAiMTk4MC0wMS0wMSIsICIxOTgwLTA3LTAxIiwKICAiMTk4MS0wNy0wMSIsICIxOTgyLTExLTAxIiwKICAiMTk5MC0wNy0wMSIsICIxOTkxLTAzLTAxIiwKICAiMjAwMS0wMy0wMSIsICIyMDAxLTExLTAxIiwKICAiMjAwNy0xMi0wMSIsICIyMDA5LTA2LTAxIiwKICAiMjAyMC0wMi0wMSIsICIyMDIwLTA1LTAxIgogICkgJT4lIAogIG11dGF0ZShQZWFrID0geW1kKFBlYWspLAogICAgICAgICBUcm91Z2ggPSB5bWQoVHJvdWdoKSkKCgpyZWNlc3Npb25fc2hhZGUgPC0gCiAgZ2VvbV9yZWN0KGRhdGEgPSByZWNlc3Npb25zLCAKICAgICAgICAgICAgaW5oZXJpdC5hZXMgPSBGLCAKICAgICAgICAgICAgYWVzKHhtaW4gPSBQZWFrLCAKICAgICAgICAgICAgICAgIHhtYXggPSBUcm91Z2gsIAogICAgICAgICAgICAgICAgeW1pbiA9IC1JbmYsIAogICAgICAgICAgICAgICAgeW1heCA9ICtJbmYpLCAKICAgICAgICAgICAgZmlsbCA9ICdwaW5rJywgCiAgICAgICAgICAgIGFscGhhID0gMC41KQoKCmhwaV9wY3QgJT4lIAogIGdncGxvdChhZXMoeCA9IHltZChkYXRlKSwgeSA9IHBjdF9jaGFuZ2UsIGNvbG9yID0gZGl2aXNpb24pKSArCiAgcmVjZXNzaW9uX3NoYWRlICsKICBnZW9tX2xpbmUoKSArCiAgdGhlbWVfbWluaW1hbCgpICsKICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDkwLCBoanVzdCA9IDEpLAogICAgICAgIHBsb3QudGl0bGUgPSBlbGVtZW50X3RleHQoaGp1c3QgPSAwLjUpLAogICAgICAgIHBsb3Quc3VidGl0bGUgPSBlbGVtZW50X3RleHQoaGp1c3QgPSAwLjUpLAogICAgICAgIHBsb3QuY2FwdGlvbiA9IGVsZW1lbnRfdGV4dChoanVzdD0wKSkgKwogIHlsYWIoIiIpICsKICB4bGFiKCJQZXJjZW50IGNoYW5nZSwgbW9udGhseSIpICsKICBnZ3RpdGxlKCJIb3VzaW5nIFByaWNlIEFwcHJlY2lhdGlvbiIsIAogICAgICAgICAgc3VidGl0bGUgPSAiYnkgVVMgQ2Vuc3VzIERpdmlzaW9uIikgKwogIGxhYnMoY2FwdGlvbiA9ICJkYXRhIHNvdXJjZTogRkhGQSIpICsKICBzY2FsZV94X2RhdGUobGltaXRzID0gYyhhcy5EYXRlKG1pbihocGlfcGN0JGRhdGUpKSwgYXMuRGF0ZShtYXgoaHBpX3BjdCRkYXRlKSkpKSAKYGBgCgoKIyBgdGltZXRrYCBmb3IgdGltZSBzZXJpZXMgZGF0YQoKVGhlIFtgdGltZXRrYF0oaHR0cHM6Ly9idXNpbmVzcy1zY2llbmNlLmdpdGh1Yi5pby90aW1ldGsvaW5kZXguaHRtbCkgcGFja2FnZSBpbmNsdWRlcyBhIGJ1bmNoIG9mIGZ1bmN0aW9ucyB0aGF0IG1ha2Ugd29ya2luZyB3aXRoIHRpbWUgc2VyaWVzIGRhdGEgc3VwZXIgZWFzeS4gVGhpcyBpbmNsdWRlcyBmdW5jdGlvbnMgZm9yIGVhc2lseSBjcmVhdGluZyBncmVhdCBsb29raW5nIHBsb3RzIG9mIHRpbWUgc2VyaWVzIGRhdGE6CgpgYGB7cn0KaHBpX3BjdCAlPiUgCiAgdW5ncm91cCgpICU+JSAKICBwbG90X3RpbWVfc2VyaWVzKGRhdGUsIHBjdF9jaGFuZ2UsIC5jb2xvcl92YXIgPSBkaXZpc2lvbiwgLnNtb290aCA9IEZBTFNFKQpgYGAKCgojIyBBbm9tYWx5IGRpYWdub3N0aWNzIHdpdGggYHRpbWV0a2AKCmB0aW1ldGtgIGFsc28gaW5jbHVkZXMgZnVuY3Rpb25zIGZvciBhdXRvbWF0aWMgYW5vbWFseSBkZXRlY3Rpb246CgpgYGB7cn0KaHBpX3BjdCAlPiUgCiAgZmlsdGVyKGRpdmlzaW9uID09ICJzb3V0aF9hdGxhbnRpY19zYSIpICU+JSAKICB1bmdyb3VwKCkgJT4lIAogIHBsb3RfYW5vbWFseV9kaWFnbm9zdGljcyhkYXRlLCBwY3RfY2hhbmdlKQpgYGAKCgojIFBsb3RseSBhbmQgYGdncGxvdGx5KClgCgpgdGltZXRrYCdzIGludGVyYWN0aXZlIHBsb3RzIHJlbHkgb24gW3Bsb3RseV0oaHR0cHM6Ly9wbG90bHkuY29tL3IvKSwgYSBsaWJyYXJ5IGZvciBidWlsZGluZyBpbnRlcmFjdGl2ZSBKYXZhU2NyaXB0IHZpc3VhbGl6YXRpb25zLiBQbG90bHkgaXMgc3VwcG9ydGVkIGluIHNldmVyYWwgZGlmZmVyZW50IGxhbmd1YWdlcyAoaW5jbHVkaW5nIFIpIGFuZCBoYXMgaXRzIG93biBzeW50YXguCgpJbXBvcnRhbnRseSwgdGhlIGBwbG90bHlgIFIgcGFja2FnZSBpbmNsdWRlcyBhIGZ1bmN0aW9uIGNhbGxlZCBgZ2dwbG90bHkoKWAgdGhhdCAoeW91IGd1ZXNzZWQgaXQhKSB0dXJucyBnZ3Bsb3RzIGludG8gcGxvdGx5IGNoYXJ0cy4KCmBgYHtyfQpocGlfcGxvdCA8LSBnZ3Bsb3QoaHBpX3RpZHksIGFlcyh4ID0gZGF0ZSwgeSA9IGhwaSwgY29sb3IgPSBkaXZpc2lvbikpICsKICBnZW9tX2xpbmUoKQoKaHBpX3Bsb3QgJT4lIAogIGdncGxvdGx5KCkKYGBgCgoKCgoKCg==